Previously, I sent around some images comparing the alpha diversity results of each run. This met with general approval as well as questions including gamma diversity, the burstiness of arrivals (co-establishment), the burstiness of extinctions (co-extinctions), and comparisons with purely neutral models of arrival and extinction.
In this document, we continue to attempt to address Chris’s question of how increased spatial mobility changes biodiversity measures. We begin by loading up each data set referring to different neutral (arrival-extinction) and spatial (dispersal speed) conditions. We can then extract our diversity metrics for each system and each site/ecosystem/environment within the system and compare them.
deSolve package does have functions that return lists, where the first element is a vector of derivatives of \(y\) with respect to \(t\). The ‘’next’’ elements seem like they could be things like the trophic structure and trophic levels.library(dplyr) # Data Manipulation
library(tidyr) # Data Pivotting
library(ggplot2) # 2-D Plot
library(plotly) # 3-D Plot
library(ggfortify) # used for biplots of PCAs
library(vegan) # Ecological analysis mega-package
library(htmltools) # Programmatic chunks
library(fitdistrplus) # Advanced fitting over MASS.
library(poweRlaw) # Fitting heavy-tailed distributions.
library(RMTRCode2) # Personal package.
# https://stackoverflow.com/a/7172832
ifrm <- function(obj, env = globalenv()) {
obj <- deparse(substitute(obj))
if(exists(obj, envir = env)) {
rm(list = obj, envir = env)
}
}
log0 <- function(x, base = exp(1)) {
xneq0 <- x != 0 & !is.na(x)
x[xneq0] <- log(x[xneq0], base = base)
return(x)
}
As of the writing of this script, all of the data is stored in the same location as this script. The parameters used to generate the interaction matrices, pools, and events is not stored with the attempts, which is something I should probably fix in the future. (I could store them in the ellipsis argument or as a separate object.)
by_for_thinning <- 10 # time steps
divide_time_by <- 1E4 # time units
burn_in <- 1E4 # time units
load_safe_thin <- function(fname, bythin, divtime, burn) {
loaded <- tryCatch({load(fname)},
error = function(e) {
print(fname)
print(e)
return(NA)
})
if (is.na(loaded)) {
return(NA)
} else {
loaded <- get(loaded)
loaded$Abundance <- loaded$Abundance[seq(from = 1,
to = nrow(loaded$Abundance),
by = bythin), ]
toEliminate <- loaded$Abundance[, -1] < loaded$Parameters$EliminationThreshold & loaded$Abundance[, -1] > 0
loaded$Abundance[, -1][toEliminate] <- 0
loaded$Abundance <- loaded$Abundance[
loaded$Abundance[, 1] > burn,
]
loaded$Abundance[, 1] <- loaded$Abundance[, 1] / divtime
return(loaded)
}
}
# All .RData
files_dat <- dir(pattern = ".RData$")
# Remove PoolMats
files_dat <- files_dat[
!grepl(x = files_dat,
pattern = "PoolMats",
fixed = TRUE)
]
# Technically overkill, but prevents unintentional loads.
# Break into two separate runs to load only intended.
# Process "MNA-FirstAttempt#####-Result-Env10-####.RData"
files_dat_FA <- files_dat[
grepl(x = files_dat,
pattern = "FirstAttempt",
fixed = TRUE)
]
Results <- sapply(
files_dat_FA,
load_safe_thin,
bythin = by_for_thinning,
divtime = divide_time_by,
burn = burn_in,
simplify = FALSE, USE.NAMES = TRUE
)
# Process "MNA-Dist#####-Ext###-Env10-####.RData"
files_dat_Ext <- files_dat[
grepl(x = files_dat,
pattern = "Dist",
fixed = TRUE)
]
Results <- c(Results, sapply(
files_dat_Ext,
load_safe_thin,
bythin = by_for_thinning,
divtime = divide_time_by,
burn = burn_in,
simplify = FALSE, USE.NAMES = TRUE
))
load_safe <- function(fname) {
loaded <- tryCatch({load(fname)},
error = function(e) {
print(fname)
print(e)
return(NA)
})
if (all(is.na(loaded))) {
return(NA)
} else {
return(sapply(loaded, get,
envir = sys.frame(sys.parent(0)),
simplify = FALSE, USE.NAMES = TRUE))
}
}
# All .RData
files_dat_PM <- dir(pattern = ".RData$")
# Remove PoolMats
files_dat_PM <- files_dat_PM[
grepl(x = files_dat_PM,
pattern = "PoolMats",
fixed = TRUE)
]
PoolsMats <- sapply(
files_dat_PM,
load_safe,
simplify = FALSE, USE.NAMES = TRUE
)
EliminiationThreshold <- unique(sapply(
Results,
function(lst) {lst$Parameters$EliminationThreshold}
))
stopifnot(length(EliminiationThreshold) == 1)
NumEnvironments <- unique(sapply(
Results,
function(lst) {lst$NumEnvironments}
))
stopifnot(length(NumEnvironments) == 1)
TrophicFunctions <- sapply(
PoolsMats,
function(PM, NE, ET) {
RMTRCode2::CalculateTrophicStructure(
Pool = PM$Pool,
NumEnvironments = NE,
InteractionMatrices = PM$InteractionMatrices,
EliminationThreshold = ET
)
},
NE = NumEnvironments,
ET = EliminiationThreshold
)
#TODO Fix warning (binding character and factor vector, coercing into character vector in bind_rows)
TrophicAnalyses <- list()
Result = Results
Nm = names(Results)
TrophFN = TrophicFunctions
for (i in seq_along(Results)) {
print(i); print(Nm[i])
# Identify Appropriate Function.
Unusual1 <- grepl(pattern = "FirstAttempt",
x = Nm[i], fixed = TRUE)
NormalKey <- strsplit(Nm[i], split = '-')[[1]][3]
linkedFN <- if(Unusual1) {
TrophFN[grepl(x = names(TrophFN),
pattern = "FirstAttempt", fixed = TRUE)][[1]]
} else {
TrophFN[grepl(x = names(TrophFN),
pattern = NormalKey, fixed = TRUE)][[1]]
}
# Apply to each row.
TrophicAnalyses[[i]] <- apply(Result[[i]]$Abundance[, -1],
MARGIN = 1,
FUN = linkedFN)
}
# TrophicAnalyses <- lapply(
# seq_along(Results),
# function(i, Result, Nm, TrophFN) {
# # Identify Appropriate Function.
# Unusual1 <- grepl(pattern = "FirstAttempt",
# x = Nm[i], fixed = TRUE)
# NormalKey <- strsplit(Nm[i], split = '-')[[1]][3]
# linkedFN <- if(Unusual1) {
# TrophFN[grepl(x = names(TrophFN),
# pattern = "FirstAttempt", fixed = TRUE)][[1]]
# } else {
# TrophFN[grepl(x = names(TrophFN),
# pattern = NormalKey, fixed = TRUE)][[1]]
# }
#
# # Apply to each row.
#
# return(apply(Result[[i]]$Abundance[, -1],
# MARGIN = 1,
# FUN = linkedFN)
# )
# },
# Result = Results,
# Nm = names(Results),
# TrophFN = TrophicFunctions
# )
names(TrophicAnalyses) <- names(Results)
# Borrowing code from FirstAttempt-Doc-Analysis.Rmd
Calculate_Diversity <- function(result) {
Diversity <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
richness <- rowSums(env != 0)
abundSum <- rowSums(env)
#NOTE: THIS CAN YIELD NAN'S (0/0).
# THIS IS NOT NECESSARILY A PROBLEM.
# IT MIGHT BE WORTH IT JUST TO USE 0 OR
# TO CATCH IT EXPLICITLY AND REPLACE WITH NAN.
entropy <- env / abundSum
entropy <- - apply(
entropy, MARGIN = 1,
FUN = function(x) {
sum(x * log0(x))
})
species <- apply(
env, MARGIN = 1,
FUN = function(x) {
toString(which(x > 0))
}
)
evenness <- entropy / log(richness)
data.frame(Time = time,
Richness = richness,
Entropy = entropy,
Evenness = evenness,
Species = species,
Environment = i,
stringsAsFactors = FALSE)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
Diversity <- dplyr::bind_rows(Diversity)
Diversity_alpha <- Diversity
# Diversity_alpha <- Diversity_alpha %>% dplyr::mutate(
# Evenness = Entropy / log(Richness)
# )
# Modify to do the gamma bits right here.
Diversity_gamma <- Diversity %>% dplyr::group_by(
Time
) %>% dplyr::summarise(
Mean = mean(Richness),
SpeciesTotal = toString(sort(unique(unlist(strsplit(paste(
Species, collapse = ", "), split = ", ", fixed = TRUE))))),
Gamma = unlist(lapply(strsplit(
SpeciesTotal, split = ", ", fixed = TRUE), function(x) length(x[x!=""]) ))
) %>% tidyr::pivot_longer(
cols = c(Mean, Gamma),
names_to = "Aggregation",
values_to = "Richness"
)
# Combine the two types of results
Diversity_alpha <- Diversity_alpha %>% dplyr::select(
-Species
) %>% tidyr::pivot_longer(
cols = c(Richness, Entropy, Evenness),
names_to = "Measurement",
values_to = "Value"
) %>% dplyr::mutate(
Environment = as.character(Environment)
)
Diversity_gamma <- Diversity_gamma %>% dplyr::select(
-SpeciesTotal
) %>% dplyr::rename(
Environment = Aggregation,
Value = Richness
) %>% dplyr::mutate(
Measurement = "Richness"
)
Diversity_beta <- Diversity_alpha %>% dplyr::filter(
Measurement == "Richness"
) %>% dplyr::select(
-Measurement
) %>% dplyr::left_join(
y = Diversity_gamma %>% dplyr::filter(
Measurement == "Richness", Environment == "Gamma"
) %>% dplyr::select(
-Measurement, -Environment
),
by = "Time",
suffix = c("_Alpha", "_Gamma")
# ) %>% dplyr::group_by(
# Time
) %>% dplyr::mutate(
BetaSpeciesMissing = Value_Gamma - Value_Alpha,
BetaSpeciesPercentage = Value_Alpha/Value_Gamma
) %>% dplyr::select(
-Value_Gamma, -Value_Alpha
) %>% tidyr::pivot_longer(
names_to = "Measurement",
values_to = "Value",
cols = c(BetaSpeciesMissing, BetaSpeciesPercentage)
# ) %>% dplyr::ungroup(
)
#print(c(colnames(Diversity_alpha), colnames(Diversity_beta), colnames(Diversity_gamma)))
Diversity <- rbind(
Diversity_alpha,
Diversity_beta,
Diversity_gamma
)
return(Diversity)
}
Diversity_jaccard_space <- function(result) {
apply(
result$Abundance,
MARGIN = 1, # Rows
function(row, envs) {
time <- row[1]
dists <- vegan::vegdist(
method = "jaccard",
x = matrix(row[-1] > 0, nrow = envs, byrow = TRUE)
)
dataf <- expand.grid(
Env1 = 1:envs,
Env2 = 1:envs
) %>% dplyr::filter(
Env1 < Env2
) %>% dplyr::mutate(
Time = time,
Jaccard = dists
)
return(dataf)
},
envs = result$NumEnvironments
)
}
Diversity_jaccard_time <- function(result, subsample = 100) {
# Break into environments, then apply it to the time series.
patches <- lapply(
1:result$NumEnvironments, function(i, abund, envs, spec) {
abund <- abund[seq(from = 1, to = nrow(abund), by = subsample), ]
times <- abund[, 1]
patch <- abund[, 1 + 1:spec + spec * (i - 1)]
dists <- vegan::vegdist(
method = "jaccard",
x = patch > 0
)
dataf <- expand.grid(
Time1 = times,
Time2 = times
) %>% dplyr::filter(
Time1 < Time2
) %>% dplyr::mutate(
Environment = i,
Jaccard = dists
)
return(dataf)
}, abund = result$Abundance, envs = result$NumEnvironments,
spec = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
}
Calculate_Species <- function(result, bintimes = FALSE) {
SpeciesPerEnvironment <- lapply(
1:result$NumEnvironments,
function(i, abund, numSpecies) {
time <- abund[, 1]
env <- abund[, 1 + 1:numSpecies + numSpecies * (i - 1)]
# Need to retrieve Position and Value
species <- apply(
cbind(time, env), MARGIN = 1,
FUN = function(x) {
time <- x[1]
dat <- x[-1]
if (any(dat > 0)) {
positions <- (which(dat > 0))
values <- dat[positions]
data.frame(
Time = time,
Species = positions,
Abundance = values,
row.names = NULL
)
} else {NULL}
# Returns as list
}
)
return(
dplyr::bind_rows(species) %>% dplyr::mutate(
Environment = i
)
)
},
abund = result$Abundance,
numSpecies = (ncol(result$Abundance) - 1) / result$NumEnvironments
)
if (bintimes) {
# Should equalise time steps.
SpeciesPerEnvironment <- lapply(
SpeciesPerEnvironment, function(SPE) {
SPE %>% dplyr::mutate(
TimeFloor = floor(Time*10)/10
) %>% dplyr::group_by(
TimeFloor, Species, Environment
) %>% dplyr::summarise(
Abundance = median(Abundance, na.rm = TRUE)
)
})
}
return(dplyr::bind_rows(SpeciesPerEnvironment))
}
# Note that if a file fails to load, we might have NA instead of a result to work with.
Diversity <- sapply(
USE.NAMES = TRUE, simplify = FALSE,
Results, function(result) {
if (length(result) == 1 && is.na(result)) {
# Problem case.
return(NA)
}
# print(paste("Calculating", Sys.time()))
# Calculate the diversity.
# We will need to extract the system properties from
# the file names which we carry through using sapply.
return(Calculate_Diversity(result))
}
)
# Expect warnings since we have all 0 rows on occasion.
JaccardSpace <- sapply(
USE.NAMES = TRUE, simplify = FALSE,
Results, function(result) {
if (length(result) == 1 && is.na(result)) {
# Problem case.
return(NA)
}
# print(paste("Calculating", Sys.time()))
# Calculate the diversity.
# We will need to extract the system properties from
# the file names which we carry through using sapply.
return(Diversity_jaccard_space(result))
}
)
# Expect warnings since we have all 0 rows on occasion.
JaccardTime <- sapply(
USE.NAMES = TRUE, simplify = FALSE,
Results, function(result) {
if (length(result) == 1 && is.na(result)) {
# Problem case.
return(NA)
}
print(paste("Calculating", Sys.time()))
# Calculate the diversity.
# We will need to extract the system properties from
# the file names which we carry through using sapply.
return(Diversity_jaccard_time(result))
# Too many time points for calculations.
# Needs to have log(length(times)^2, base = 2) < 31.
# => length(times) < 46340 or so.
}
)
[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:55"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:56"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:57"
[1] "Calculating 2021-11-14 18:47:58"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:47:59"
[1] "Calculating 2021-11-14 18:48:00"
[1] "Calculating 2021-11-14 18:48:00"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:01"
[1] "Calculating 2021-11-14 18:48:02"
[1] "Calculating 2021-11-14 18:48:02"
SpeciesPresence <- sapply(
USE.NAMES = TRUE, simplify = FALSE,
Results, function(result) {
if (length(result) == 1 && is.na(result)) {
# Problem case.
return(NA)
}
# print(paste("Calculating", Sys.time()))
# Calculate the diversity.
# We will need to extract the system properties from
# the file names which we carry through using sapply.
return(Calculate_Species(result))
}
)
Properties <- strsplit(names(Diversity), '-',
fixed = TRUE)
# 1st Chunk: Name, Discard
# 2nd Chunk: Iteration + Distance
# 3rd Chunk: Result or Extinction Rate (or Arrival?)
# 4th Chunk: Number of Environments
# 5th Chunk: Space Type + .RData
# Note the mix of Keyword and Location Structure (oops).
# Note also that this strsplit character is a bad decision and should be changed for next time. (D'oh.)
# (E.g. Dates DD-MM-YYYY, Decimals 1.35e-05.)
Properties <- data.frame(
do.call(rbind, Properties),
stringsAsFactors = FALSE
)
names(Properties)[1:5] <- c(
"Name", "IterANDDist", "Modifier", "EnvNum", "SpaceAND.RData"
)
Properties$FullName <- names(Diversity)
# Capture the position between the text (first group)
# and the set of numbers (somehow without the +).
# The \\K resets so that we do not capture any text.
patternString <- "((?>[a-zA-Z]+)(?=[0-9eE]))\\K"
# Split strings. Some of the trick will be to introduce
# a character to make the separation around. We use "_".
Properties <- Properties %>% dplyr::mutate(
IterANDDist = gsub(pattern = patternString,
replacement = "_",
x = IterANDDist, perl = TRUE),
Modifier = gsub(pattern = patternString,
replacement = "_",
x = Modifier, perl = TRUE),
EnvNum = gsub(pattern = patternString,
replacement = "_",
x = EnvNum, perl = TRUE)
) %>% tidyr::separate(
IterANDDist, into = c("Iter", "Distance"),
sep = "[_]", fill = "right"
) %>% tidyr::separate(
Modifier, into = c("Modifier", "ModIntensity"),
sep = "[_]", fill = "right"
) %>% tidyr::separate(
EnvNum, into = c("Env", "Environments"),
sep = "[_]"
) %>% tidyr::separate(
SpaceAND.RData, into = c("Space", ".RData"),
sep = "[.]"
) %>% dplyr::select(
-Name, -.RData, -Env
) %>% dplyr::mutate(
Distance = dplyr::case_when(
is.na(Distance) ~ "1e+00",
TRUE ~ Distance
)
)
Diversity <- lapply(1:length(Diversity),
function(i, df, nm) {
df[[i]] %>% mutate(
Simulation = nm[i]
)
},
df = Diversity,
nm = names(Diversity))
JaccardSpace <- lapply(1:length(JaccardSpace),
function(i, df, nm) {
df[[i]] %>% dplyr::bind_rows(
# Need to account for the by times...
# Generates attribute warnings.
) %>% mutate(
Simulation = nm[i]
)
},
df = JaccardSpace,
nm = names(JaccardSpace))
JaccardTime <- lapply(1:length(JaccardTime),
function(i, df, nm) {
df[[i]] %>% dplyr::bind_rows(
# Need to account for the by envs...
# Generates attribute warnings.
) %>% mutate(
Simulation = nm[i]
)
},
df = JaccardTime,
nm = names(JaccardTime))
SpeciesPresence <- lapply(1:length(SpeciesPresence),
function(i, df, nm) {
df[[i]] %>% mutate(
Simulation = nm[i]
)
},
df = SpeciesPresence,
nm = names(SpeciesPresence))
Diversity <- dplyr::left_join(
dplyr::bind_rows(Diversity),
Properties,
by = c("Simulation" = "FullName")
)
JaccardSpace <- dplyr::left_join(
dplyr::bind_rows(JaccardSpace),
Properties,
by = c("Simulation" = "FullName")
)
JaccardTime <- dplyr::left_join(
dplyr::bind_rows(JaccardTime),
Properties,
by = c("Simulation" = "FullName")
)
SpeciesPresence <- dplyr::left_join(
dplyr::bind_rows(SpeciesPresence),
Properties,
by = c("Simulation" = "FullName")
)
Diversity <- Diversity %>% dplyr::mutate(
Distance = dplyr::case_when(
is.na(Distance) ~ "1e+00",
TRUE ~ Distance
)
)
ggplot2::ggplot(
Diversity %>% dplyr::filter(
Measurement == "Richness",
Environment != "Gamma",
Space != "Ring",
Environment != "Mean"
),
ggplot2::aes(
x = Time,
y = Value,
color = factor(Environment),
alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::filter(
Measurement == "Richness",
Environment != "Gamma",
Space != "Ring",
Environment == "Mean"
),
color = "black"
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
) + ggplot2::labs(
y = "Richness",
x = paste0("Time, ", divide_time_by, " units"),
title = "Overview of Alpha Richness over Time by System Properties",
caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
)
# ggplot2::ggsave(overallrich + ggplot2::coord_cartesian(ylim = c(0, 25)), filename = "MNA-AlphaRichness-Overview.pdf", dpi = "retina", width = 11, height = 8)
Plot_Richness <- function(df) {
tempname <- paste(unique(df$Simulation), collapse = " ")
temp <- ggplot2::ggplot(
df %>% dplyr::filter(
Measurement == "Richness",
Environment != "Gamma",
Environment != "Mean"
),
ggplot2::aes(
x = Time,
y = Value,
color = factor(Environment),
alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = df %>% dplyr::filter(
Measurement == "Richness",
Environment != "Gamma",
Environment == "Mean"
),
color = "black"
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
) + ggplot2::labs(
title = tempname,
x = paste0("Time, ", divide_time_by, " units"),
y = "Richness"
)
return(temp)
}
Plots_Richness_alpha <- Diversity %>% dplyr::group_split(
Simulation
) %>% purrr::map(
Plot_Richness
)
# In the next chunk, we use
# www.r-bloggers.com/2020/07/programmatically-create-new-geadings-and-outputs-in-rmarkdown/
ggplot2::ggplot(
Diversity %>% dplyr::filter(
Measurement == "Evenness",
Environment != "Gamma",
Space != "Ring",
Environment != "Mean"
),
ggplot2::aes(
x = Time,
y = Value,
color = factor(Environment),
alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::filter(
Measurement == "Evenness",
Environment != "Gamma",
Space != "Ring",
Environment == "Mean"
),
color = "black"
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
) + ggplot2::labs(
y = "Evenness",
x = paste0("Time, ", divide_time_by, " units"),
title = "Overview of Evenness over Time by System Properties",
caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
)
# ggplot2::ggsave(overallrich + ggplot2::coord_cartesian(ylim = c(0, 25)), filename = "MNA-AlphaRichness-Overview.pdf", dpi = "retina", width = 11, height = 8)
Plot_Evenness <- function(df) {
tempname <- paste(unique(df$Simulation), collapse = " ")
temp <- ggplot2::ggplot(
df %>% dplyr::filter(
Measurement == "Evenness",
Environment != "Gamma",
Environment != "Mean"
),
ggplot2::aes(
x = Time,
y = Value,
color = factor(Environment),
alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = df %>% dplyr::filter(
Measurement == "Evenness",
Environment != "Gamma",
Environment == "Mean"
),
color = "black"
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
) + ggplot2::labs(
title = tempname,
x = paste0("Time, ", divide_time_by, " units"),
y = "Evenness"
)
return(temp)
}
Plots_Evenness <- Diversity %>% dplyr::group_split(
Simulation
) %>% purrr::map(
Plot_Evenness
)
# In the next chunk, we use
# www.r-bloggers.com/2020/07/programmatically-create-new-headings-and-outputs-in-rmarkdown/
ggplot2::ggplot(
Diversity %>% dplyr::filter(
Measurement == "Richness",
Environment == "Gamma",
Space != "Ring"
),
ggplot2::aes(
x = Time,
y = Value
)
) + ggplot2::geom_line(
) + ggplot2::labs(
y = "Richness",
x = paste0("Time, ", divide_time_by, " units"),
title = "Gamma Richness over Time by System Properties",
caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::coord_cartesian(ylim = c(0, 45))
# ggplot2::ggsave(overallgamma, filename = "MNA-GammaRichness-Overview.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
Diversity %>% dplyr::filter(
Measurement == "BetaSpeciesMissing",
Environment != "Gamma",
Space != "Ring",
Environment != "Mean"
),
ggplot2::aes(
x = Time,
y = Value,
color = factor(Environment),
alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::filter(
Measurement == "BetaSpeciesMissing",
Environment != "Gamma",
Space != "Ring",
Environment != "Mean"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance, Time
) %>% dplyr::summarise(
Environment = "Mean",
Value = mean(Value, na.rm = TRUE)
), color = "black"
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
) + ggplot2::labs(
y = "Absolute Species Turnover", # en.wikipedia.org/wiki/Beta_diversity
x = paste0("Time, ", divide_time_by, " units"),
title = "Overview of Gamma - Alpha Richness over Time by System Properties",
caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
)
# ggplot2::ggsave(overallbetamiss, filename = "MNA-BetaMissing-Overview.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
Diversity %>% dplyr::filter(
Measurement == "BetaSpeciesPercentage",
Environment != "Gamma",
Space != "Ring",
Environment != "Mean"
),
ggplot2::aes(
x = Time,
y = Value,
color = factor(Environment),
alpha = ifelse(Environment %in% c("3", "7", "Mean"), 1, 0.3)
)
) + ggplot2::geom_line(
) + ggplot2::geom_line(
data = Diversity %>% dplyr::filter(
Measurement == "BetaSpeciesPercentage",
Environment != "Gamma",
Space != "Ring",
Environment != "Mean"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance, Time
) %>% dplyr::summarise(
Environment = "Mean",
Value = mean(Value, na.rm = TRUE)
), color = "black"
) + ggplot2::guides(
alpha = "none"
) + ggplot2::scale_color_discrete(
"Environment"
) + ggplot2::labs(
y = "Percentage Species Present", # en.wikipedia.org/wiki/Beta_diversity
x = paste0("Time, ", divide_time_by, " units"),
title = "Overview of Alpha/Gamma Richness over Time by System Properties",
caption = "Each row has the same arrival and extinction events."
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
)
# ggplot2::ggsave(overallbetapercent, filename = "MNA-BetaPercent-Overview.pdf", dpi = "retina", width = 11, height = 8)
Diversity_AG <- dplyr::left_join(
Diversity %>% dplyr::filter(
Measurement == "Richness",
Space != "Ring",
Environment != "Gamma"
),
Diversity %>% dplyr::filter(
Measurement == "Richness",
Space != "Ring",
Environment == "Gamma"
) %>% dplyr::rename(
Gamma = Value
) %>% dplyr::select(
-Environment
)
) %>% dplyr::mutate(
Distance = ifelse(is.na(Distance), 1, Distance)
)
Joining, by = c("Time", "Measurement", "Simulation", "Iter", "Distance", "Modifier", "ModIntensity", "Environments", "Space")
Diversity_AG_Binned <- Diversity_AG %>% dplyr::mutate(
TimeFloor = floor(Time * 10) / 10
) %>% dplyr::distinct(
Modifier, ModIntensity, Space, Distance, # Simulation/Facets
Environment, TimeFloor, # Grouping Bins
Value, Gamma # Values. If either moves, the trajectory entered a new square.
)
# ggplot2::ggplot(
# temp %>% dplyr::filter(
# Modifier == "Result",
# Space == "None",
# Environment != "Mean"
# ), ggplot2::aes(
# x = Gamma,
# y = Value
# )
# ) + ggplot2::geom_bin2d(
# ) + ggplot2::scale_fill_viridis_c(
# ) + ggplot2::labs(
# Title = "No Spatial Structure, Equal Extinction & Arrival Rates",
# x = "Gamma Richness",
# y = "Alpha Richness"
# )
ggplot2::ggplot(
Diversity_AG %>% dplyr::filter(
Environment != "Mean",
Space != "Ring"
), ggplot2::aes(
x = Value,
y = Gamma
)
) + ggplot2::geom_bin2d(
) + ggplot2::scale_fill_viridis_c(
trans = "log10"
) + ggplot2::geom_point(
data = Diversity_AG %>% dplyr::filter(
Environment != "Mean",
Space != "Ring"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance
) %>% dplyr::summarise(
Gamma = mean(Gamma, na.rm = TRUE),
Value = mean(Value, na.rm = TRUE)
), color = "red", size = 2, shape = 4
) + ggplot2::labs(
x = "Alpha Richness",
y = "Gamma Richness"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_abline(slope = 1, intercept = 0)
# ggplot2::ggsave(overallalphagamma, filename = "MNA-AlphaGamma-Overview.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
Diversity_AG_Binned %>% dplyr::filter(
Environment != "Mean",
Space != "Ring"
), ggplot2::aes(
x = Value,
y = Gamma
)
) + ggplot2::geom_bin2d(
) + ggplot2::scale_fill_viridis_c(
trans = "log10"
) + ggplot2::geom_point(
data = Diversity_AG_Binned %>% dplyr::filter(
Environment != "Mean",
Space != "Ring"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance
) %>% dplyr::summarise(
Gamma = mean(Gamma, na.rm = TRUE),
Value = mean(Value, na.rm = TRUE)
), color = "red", size = 2, shape = 4
) + ggplot2::labs(
x = "Alpha Richness",
y = "Gamma Richness",
subtitle = "Trajectories binned by time to equalise time steps."
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_abline(slope = 1, intercept = 0)
# ggplot2::ggsave(overallalphagammabin, filename = "MNA-AlphaGammaBinned-Overview.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
Diversity_AG_Binned %>% dplyr::filter(
Environment != "Mean",
Space != "Ring"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance
) %>% dplyr::summarise(
Gamma = mean(Gamma, na.rm = TRUE),
Value = mean(Value, na.rm = TRUE)
),
ggplot2::aes(
x = Value,
y = Gamma,
color = interaction(Modifier, ModIntensity),
shape = interaction(Space, Distance)
)
) + ggplot2::geom_point(
size = 2
) + ggplot2::labs(
x = "Alpha Richness",
y = "Gamma Richness",
subtitle = "Trajectories binned by time to equalise time steps."
) + ggplot2::geom_abline(slope = 1, intercept = 0)
# ggplot2::ggsave(overallalphagammamean, filename = "MNA-AlphaGammaBinned-Means.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
JaccardSpace %>% dplyr::filter(Time > 3, Space != "Ring"),
ggplot2::aes(x = Time, y = Jaccard,
color = interaction(Env1, Env2)
)
) + ggplot2::geom_line(
alpha = 0.3
) + ggplot2::geom_line(
data = JaccardSpace %>% dplyr::filter(Time > 3, Space != "Ring") %>% dplyr::group_by(
Time, Distance, Modifier, ModIntensity, Space
) %>% dplyr::summarise(
Jaccard = mean(Jaccard, na.rm = TRUE)
),
ggplot2::aes(
x = Time, y = Jaccard
),
inherit.aes = FALSE, color = "black"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::guides(
color = "none"
)
ggplot2::ggplot(
JaccardTime %>% dplyr::filter(
Space != "Ring"
) %>% dplyr::group_by(
Environment, Distance, Space, Modifier, ModIntensity
) %>% dplyr::mutate(
TimeDifference = Time2 - Time1
),
ggplot2::aes(
x = TimeDifference,
y = Jaccard
)
) + ggplot2::geom_bin2d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_fill_viridis_c(
trans = "log"
)
ggplot2::ggplot(
JaccardTime %>% dplyr::filter(
Space != "Ring",
Space == "None",
Modifier == "Result"
) %>% dplyr::group_by(
Environment, Distance, Space, Modifier, ModIntensity
) %>% dplyr::mutate(
TimeDifference = Time2 - Time1
),
ggplot2::aes(
x = TimeDifference,
y = Jaccard
)
) + ggplot2::geom_bin2d(
) + ggplot2::facet_wrap(
~ Environment
) + ggplot2::scale_fill_viridis_c(
trans = "log"
)
# ggplot2::ggplot(
# SpeciesPresence %>% dplyr::filter(
# Space != "Ring",
# TimeFloor > 0.2 # Remove "burn-in" which has "impossibly" high presence.
# ),
# ggplot2::aes(
# x = TimeFloor,
# y = Species
# )
# ) + ggplot2::geom_bin2d(
# binwidth = c(0.1, 1)
# ) + ggplot2::scale_fill_viridis_c(
# )+ ggplot2::facet_grid(
# Modifier + ModIntensity ~ Space + Distance
# )
ggplot2::ggplot(
SpeciesPresence %>% dplyr::filter(
Space != "Ring"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance,
Species, Time
) %>% dplyr::summarise(
Count = n()
),
ggplot2::aes(x = Time, y = Species, color = Count)
) + ggplot2::geom_point(
shape = '.'
) + ggplot2::scale_color_viridis_c(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
yintercept = 34.5, color = "red"
)
ggplot2::ggplot(
SpeciesPresence %>% dplyr::filter(
Space != "Ring"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance,
Species, Time
) %>% dplyr::summarise(
Count = n()
),
ggplot2::aes(x = Time, y = Species, color = Count)
) + ggplot2::geom_point(
shape = '.'
) + ggplot2::scale_color_viridis_c(
direction = -1
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
yintercept = 34.5, color = "red"
)
# Brute forcing here because I would rather position myself to move on if possible.
# This works because all of the pools are the same (except that the original does not use factors because stringsAsFactors differed between the machines generating the systems).
# > all.equal(PoolsMats$`MNA-Ext10-PoolMats-Env10.RData`, PoolsMats$`MNA-Ext0.1-PoolMats-Env10.RData`)
# [1] TRUE
# > all.equal(PoolsMats$`MNA-Ext10-PoolMats-Env10.RData`, PoolsMats$`MNA-Arr0.1-PoolMats-Env10.RData`)
# [1] TRUE
# > all.equal(PoolsMats$`MNA-Ext10-PoolMats-Env10.RData`, PoolsMats$`MNA-Arr10-PoolMats-Env10.RData`)
# [1] TRUE
SpeciesPresence$Sizes <- PoolsMats$`MNA-FirstAttempt-PoolMats-Env10.RData`$Pool$Size[SpeciesPresence$Species]
ggplot2::ggplot(
SpeciesPresence %>% dplyr::filter(
Space != "Ring"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance,
Species, Time, Sizes
) %>% dplyr::summarise(
Count = n()
) %>% dplyr::arrange(
Sizes
) %>% dplyr::ungroup(
) %>% dplyr::mutate(
Species = factor(Species, levels = Species)
),
ggplot2::aes(x = Time, y = Species, color = Count)
) + ggplot2::geom_point(
shape = '.'
) + ggplot2::scale_color_viridis_c(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
yintercept = 34.5, color = "red"
)
Error in `levels<-`(`*tmp*`, value = as.character(levels)) :
factor level [2] is duplicated
# ggplot2::ggplot(
# SpeciesPresence %>% dplyr::filter(
# Space != "Ring",
# TimeFloor > 0.2 # Remove "burn-in" which has "impossibly" high presence.
# ),
# ggplot2::aes(
# x = TimeFloor,
# y = Abundance
# )
# ) + ggplot2::geom_bin2d(
# ) + ggplot2::scale_fill_viridis_c(
# ) + ggplot2::facet_grid(
# Modifier + ModIntensity ~ Space + Distance
# ) + ggplot2::scale_y_log10(
# )
ggplot2::ggplot(
SpeciesPresence %>% dplyr::filter(
Space != "Ring"
) %>% dplyr::mutate(
Abundance = round(Abundance)
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance,
Abundance, Time
) %>% dplyr::summarise(
Count = n()
),
ggplot2::aes(x = Time, y = Abundance, color = Count)
) + ggplot2::geom_point(
shape = '.'
) + ggplot2::scale_color_viridis_c(
trans = "log"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_y_log10(
)
(Recall that there is a threshold above which we have consumers and below which we have producers. This will be denoted with a bright line at the threshold.)
# ggplot2::ggplot(
# SpeciesPresence %>% dplyr::filter(
# Space != "Ring",
# TimeFloor > 0.2 # Remove "burn-in" which has "impossibly" high presence.
# ),
# ggplot2::aes(
# x = TimeFloor,
# y = Sizes
# )
# ) + ggplot2::geom_bin2d(
# binwidth = c(0.1, 0.03)
# ) + ggplot2::scale_fill_viridis_c(
# ) + ggplot2::facet_grid(
# Modifier + ModIntensity ~ Space + Distance
# ) + ggplot2::scale_y_log10(
# ) + ggplot2::geom_hline(
# yintercept = 0.1, color = "red"
# )
ggplot2::ggplot(
SpeciesPresence %>% dplyr::filter(
Space != "Ring"
) %>% dplyr::group_by(
Modifier, ModIntensity, Space, Distance,
Sizes, Time
) %>% dplyr::summarise(
Count = n()
),
ggplot2::aes(x = Time, y = Sizes, color = Count)
) + ggplot2::geom_point(
shape = '.'
) + ggplot2::scale_color_viridis_c(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::geom_hline(
yintercept = 0.1, color = "red"
) + ggplot2::scale_y_log10(
)
Compare with the distribution of traits amongst species.
ggplot2::ggplot(
PoolsMats$`MNA-FirstAttempt-PoolMats-Env10.RData`$Pool,
ggplot2::aes(
x = Size,
fill = Type
)
) + ggplot2::geom_histogram(#ggplot2::geom_density(
#alpha = 0.5, adjust = 1/10
bins = 100
) + ggplot2::scale_x_log10() + ggplot2::coord_flip()
Calculate_Events <- function(result, abundanceTime = divide_time_by, thinning = by_for_thinning) {
abund <- result$Abundance
abund[, -1] <- abund[, -1] > 0 # Presence-Absence
abundDiff <- apply(abund[, -1], 2, diff)
# Arrival (1) Elimination (-1)
# if arrival at time 93,
# 92nd entry is 1.
abundDiff <- cbind(
abund[-1, 1] * abundanceTime, abundDiff
)
Events <- lapply(
1:ncol(abundDiff[, -1]),
function(i, ab, tm, sp) {
arrivals <- which(ab[, i] == 1)
extincts <- which(ab[, i] == -1)
if (length(arrivals) + length(extincts) == 0) {
return(NULL)
}
data.frame(
Times = c(tm[arrivals], tm[extincts]),
Species = ((i - 1) %% sp) + 1, # 1:1000 -> 1:100
Environment = ((i - 1) %/% sp) + 1,
Type = c(rep("Arrival", length(arrivals)),
rep("Extinct", length(extincts))),
stringsAsFactors = FALSE
)
},
ab = abundDiff[, -1],
tm = abundDiff[, 1],
sp = ncol(abundDiff[, -1]) / result$NumEnvironments#,
#ne = result$NumEnvironments
) %>% dplyr::bind_rows(
) %>% dplyr::arrange(
Times, Species, Environment, Type
)
# Now we check to see which events are in the record.
# Note that, due to thinning, we do have a theoretical
# problem: an event time might be misrecorded when
# extracted from the abundance record.
# Hence we cannot just use filtering join operations.
# Since we know the maximum time step size and the
# thinning, we know that we should detect a change
# within (thinning) * (max. time step size) units.
maximumGap <- thinning * result$Parameters$MaximumTimeStep
result$Events$Type <- as.character(result$Events$Type)
# Connect all same event, even with different times.
# Remove those that cannot be the same event.
# Treat this as the list of Neutral Events that
# actually happened.
EventsOfficial <- result$Events %>% dplyr::left_join(
Events,
by = c("Species", "Environment", "Type")
) %>% dplyr::mutate(
# dplyr::filter( # Filtering does not work.
# We are seeing losses of about 50% with filter.
# "Rows in x with no match in y will have NA values in the new columns."
# -> Times.y will be NA.
# Times.y is the detected time.
# Times.x is the recorded action's time.
# Note that we can get false readings from subtracting two almost the same
# numbers, so we need to appeal to machine precision.
# See all.equal's tolerance argument.
`Times.y` = dplyr::case_when(
is.na(`Times.y`) ~ as.double(NA),
(
`Times.y` - `Times.x` < maximumGap + sqrt(.Machine$double.eps) &
`Times.y` - `Times.x` >= -sqrt(.Machine$double.eps)) ~ `Times.y`,
TRUE ~ as.double(NA)
)
) %>% dplyr::group_by(
`Times.x`, Species, Environment, Type, Success
# Don't discard Success, others are true groups.
# We want to preserve the first `Times.y`
# (in case an event happens multiple times)
# but if there are no numerics,
# we instead want to keep one of the NAs.
) %>% dplyr::summarise(
`Times.y` = if (length(na.omit(`Times.y`)) == 0) NA else min(`Times.y`, na.rm = TRUE)
) %>% dplyr::ungroup(
) %>% dplyr::rename(
TimeImplemented = `Times.x`,
TimeDetected = `Times.y`
)
# Events that were not detected but were successful.
# If this is the case, something happened and was
# undone in the same timespan.
# This might happen due to arrivals from adjacent
# patches.
# I.e. elimination coinciding with arrival or
# arrival being too dissipated by dispersal.
# This probably shouldn't happen in the disconnected system.
EventsNotDetected <- EventsOfficial %>% dplyr::filter(
is.na(TimeDetected), Success == TRUE
# ) %>% dplyr::select(
# -`Times.y`
# ) %>% dplyr::rename(
# Times = `Times.x`
) %>% dplyr::mutate(
Neutral = TRUE,
Detected = FALSE
)
# Events that were detected and were successful. 'True Positives'
# We recorded them as having happened in Events and Abundance.
# These are also Neutral with high probability.
EventsDetected <- EventsOfficial %>% dplyr::filter(
!is.na(TimeDetected), Success == TRUE
# ) %>% dplyr::select(
# # We keep Times.y for comparison with Events.
# -`Times.x`
# ) %>% dplyr::rename(
# Times = `Times.y`
) %>% dplyr::mutate(
Neutral = TRUE,
Detected = TRUE
)
# Events that were not detected and were not successful. 'True Negatives'
# Also Events that were detected and were not successful. 'False Positives'?
# The detection must be of a different event if the event
# we recorded was unsuccessful after all.
EventsFailed <- EventsOfficial %>% dplyr::filter(
#is.na(`Times.y`),
Success != TRUE
# ) %>% dplyr::select(
# -`Times.y`
# ) %>% dplyr::rename(
# Times = `Times.x`
) %>% dplyr::mutate(
Neutral = TRUE,
Detected = FALSE
)
# So we have events that were detected but not successful.
# Such events should probably be listed twice:
# once as neutral (the failure, above)
# once as non-neutral (the detected event, below).
# (Note Events are from abundance and thus detected.)
EventsNotOfficial <- Events %>% dplyr::rename(
TimeDetected = Times
) %>% dplyr::anti_join(
EventsDetected, by = c("TimeDetected", "Species", "Environment", "Type")
) %>% dplyr::mutate(
Success = TRUE,
Neutral = FALSE,
Detected = TRUE
)
# The remainder of event space is events that are
# not neutral and not successful or
# events that were not implemented.
# Note then that, if everything went well
stopifnot(nrow(EventsNotOfficial) + nrow(EventsDetected) == nrow(Events),
nrow(EventsDetected) + nrow(EventsFailed) + nrow(EventsNotDetected) == nrow(result$Events))
return(dplyr::bind_rows(
EventsDetected,
EventsNotDetected,
EventsFailed,
EventsNotOfficial
) %>% dplyr::mutate(
Times = dplyr::case_when(
!is.na(TimeImplemented) ~ TimeImplemented,
!is.na(TimeDetected) ~ TimeDetected
)
) %>% dplyr::arrange(Times, Environment, Species, Type))
}
Events <- sapply(
USE.NAMES = TRUE, simplify = FALSE,
Results, function(result) {
if (length(result) == 1 && is.na(result)) {
# Problem case.
return(NA)
}
# print(paste("Calculating", Sys.time()))
# Calculate the diversity.
# We will need to extract the system properties from
# the file names which we carry through using sapply.
return(Calculate_Events(result))
}
)
Events <- lapply(1:length(Events),
function(i, df, nm) {
df[[i]] %>% mutate(
Simulation = nm[i]
)
},
df = Events,
nm = names(Events))
EventsSuccesses <- lapply(Events, function(Event) {
Event %>% dplyr::filter(
Success == TRUE
) %>% dplyr::group_by(
Environment
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
)
})
Events <- lapply(Events, function(Event) {
Event %>% dplyr::group_by(
Environment
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
)
})
Events <- dplyr::left_join(
dplyr::bind_rows(Events),
Properties,
by = c("Simulation" = "FullName")
)
EventsSuccesses <- dplyr::left_join(
dplyr::bind_rows(EventsSuccesses),
Properties,
by = c("Simulation" = "FullName")
)
Events <- Events %>% dplyr::mutate(
Distance = dplyr::case_when(
is.na(Distance) ~ "1e+00",
TRUE ~ Distance
)
)
EventsSuccesses <- EventsSuccesses %>% dplyr::mutate(
Distance = dplyr::case_when(
is.na(Distance) ~ "1e+00",
TRUE ~ Distance
)
)
# Retrieve the Characteristic Rate used.
# Since all of the Pools and Matrices are the same, they all have the same characteristic rate.
# There is a fudge here. Looking back at the code, the first matrix was used, rather than looking over all of the matrices, but we note that the results are essentially the same.
# Perhaps something to be more careful of for the next round of simulations...
CharacteristicRate <- max(abs(eigen(PoolsMats[[1]]$InteractionMatrices$Mats[[1]])$values))
set.seed(1)
exampleData <- Events[[2]] %>% dplyr::filter(Neutral) %>% dplyr::arrange(Times) %>% pull(Times) %>% diff
fitdistrplus::descdist(exampleData, boot = 100)
censData <- ifelse(exampleData == 0, .Machine$double.eps, exampleData)
exampleFitsNotHeavy <- list(
exp = fitdistrplus::fitdist(censData, "exp", method = "mle"),
gamma = fitdistrplus::fitdist(censData,
"gamma", method = "mle"),
weibull = fitdistrplus::fitdist(censData,
"weibull", method = "mle")
)
exampleFitsHeavy <- lapply(
c(poweRlaw::conexp, poweRlaw::conlnorm, poweRlaw::conpl, poweRlaw::conweibull), function(f) {f$new(censData)}
)
exampleFitsHeavy <- lapply(
exampleFitsHeavy,
function(d) {
d$setXmin(poweRlaw::estimate_xmin(d, xmax = Inf))
d
}
)
names(exampleFitsHeavy) <- c("exp", "lnorm", "pl", "weibull")
fitdistrplus::gofstat(exampleFitsNotHeavy)
quiet(exampleData_gamlss <- gamlss::fitDist(exampleData))
exampleData_gamlss
quiet(censData_gamlss <- gamlss::fitDist(censData))
exampleData_gamlss
exampleFitsHeavy
ggplot2::ggplot(
Events %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() #-> environmenteventssequence
# ggplot2::ggsave(environmenteventssequence, filename = "MNA-EventEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
Events %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() #-> environmenteventsneutral
# ggplot2::ggsave(environmenteventsneutral, filename = "MNA-EventEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)
Events %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() + ggplot2::coord_cartesian(ylim = c(0, 2)) # -> overalleventssequence
# ggplot2::ggsave(overalleventssequence, filename = "MNA-EventOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
Events %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() # -> overalleventsneutral
# ggplot2::ggsave(overalleventsneutral, filename = "MNA-EventOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
EventsSuccesses %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
ylim = c(0, 2)
) #-> environmenteventssuccesssequence
# ggplot2::ggsave(environmenteventssuccesssequence, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
EventsSuccesses %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() #-> environmenteventssuccessneutral
# ggplot2::ggsave(environmenteventssuccessneutral, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)
EventsSuccesses %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() + ggplot2::coord_cartesian(ylim = c(0, 2)) # -> overalleventssuccesssequence
# ggplot2::ggsave(overalleventssuccesssequence, filename = "MNA-EventSuccessOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
EventsSuccesses %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10() # -> overalleventssuccessneutral
# ggplot2::ggsave(overalleventssuccessneutral, filename = "MNA-EventSuccessOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
Events %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4)
) #-> environmenteventssequence
# ggplot2::ggsave(environmenteventssequence, filename = "MNA-EventEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
Events %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4)
) #-> environmenteventsneutral
# ggplot2::ggsave(environmenteventsneutral, filename = "MNA-EventEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)
Events %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4),
ylim = c(0, 2)
) # -> overalleventssequence
# ggplot2::ggsave(overalleventssequence, filename = "MNA-EventOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
Events %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4)
) # -> overalleventsneutral
# ggplot2::ggsave(overalleventsneutral, filename = "MNA-EventOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
EventsSuccesses %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4),
ylim = c(0, 2)
) #-> environmenteventssuccesssequence
# ggplot2::ggsave(environmenteventssuccesssequence, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
ggplot2::ggplot(
EventsSuccesses %>% dplyr::filter(
Space != "Ring"
),
ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4)
) #-> environmenteventssuccessneutral
# ggplot2::ggsave(environmenteventssuccessneutral, filename = "MNA-EventSuccessEnvironmentArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)
EventsSuccesses %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = factor(EventInSequence)
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
"Event #"
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4),
ylim = c(0, 2)
) # -> overalleventssuccesssequence
# ggplot2::ggsave(overalleventssuccesssequence, filename = "MNA-EventSuccessOverallArrivalTimes-Sequence.pdf", dpi = "retina", width = 11, height = 8)
EventsSuccesses %>% dplyr::group_by(
Simulation, Iter, Distance,
Modifier, ModIntensity, Environments, Space
) %>% dplyr::arrange(
Times
) %>% dplyr::mutate(
InterarrivalTime = Times - lag(Times),
Dynamic = !Neutral,
Sequence = cumsum(Neutral)
) %>% dplyr::group_by(
Environment, Sequence
) %>% dplyr::mutate(
EventInSequence = cumsum(Dynamic)
) %>% dplyr::filter(
Space != "Ring"
) %>% ggplot2::ggplot(ggplot2::aes(
x = InterarrivalTime,
fill = Neutral
)
) + ggplot2::geom_density(
alpha = 0.25
) + ggplot2::scale_fill_viridis_d(
) + ggplot2::facet_grid(
Modifier + ModIntensity ~ Space + Distance
) + ggplot2::scale_x_log10(
) + ggplot2::coord_cartesian(
xlim = c(1e-1, 1e4)
) # -> overalleventssuccessneutral
# ggplot2::ggsave(overalleventssuccessneutral, filename = "MNA-EventSuccessOverallArrivalTimes-Overview.pdf", dpi = "retina", width = 11, height = 8)